
%function impliedCalvoSlope 
clear all; 

dim.draws=50000; 

flags.Kimball=0; 
fixed.phi=1.5; 


par.betin.dist='G'; 
par.betin.mean=0.17; 
par.betin.std=0.05;
par.betin.ubnd=2; 

par.sigmac.dist='N'; 
par.sigmac.mean=1; 
par.sigmac.std=0.5;
par.sigmac.lbnd=1e-5; 
par.sigmac.ubnd=10; 

par.zss.dist='N'; 
par.zss.mean=0.5; 
par.zss.std=0.05; 
par.zss.lbnd=1e-5; 
par.zss.ubnd=1; 

par.zeta.dist='B'; 
par.zeta.mean=0.85; 
par.zeta.std=0.025; 

par.iota.dist='B'; 
par.iota.mean=0.5; 
par.iota.std=0.15; 

names.par=fieldnames(par); 
dim.par=length(names.par); 

prior.dist=cell(dim.par,1); 
prior.mean=zeros(dim.par,1); 
prior.std =zeros(dim.par,1);
prior.alpha=zeros(dim.par,1); 
prior.beta =zeros(dim.par,1); 

for ii=1:dim.par;
    fldnm=char(names.par{ii});
    prior.dist(ii)={par.(fldnm).dist};
    prior.mean(ii)=par.(fldnm).mean;
    prior.std(ii)=par.(fldnm).std;
    
    if isempty(strmatch(char(prior.dist(ii)),{'G','B','I','W'}))==false; 
        prior.lbnd(ii)=1e-5; 
    else 
        prior.lbnd(ii)=par.(fldnm).lbnd;
    end
    
    if strcmpi( char(prior.dist(ii)) , 'B' ) == true
        prior.ubnd(ii)=1;
    else
        prior.ubnd(ii)=par.(fldnm).ubnd;
    end
            
    pos.(fldnm)=ii; 
end

% prior.mean(ii)=par.(fldnm).
% prior.mean(ii)=par.(fldnm).    

[prior.alpha,prior.beta]=msdtoab(prior.mean,prior.std,prior.dist); 
drawsmat=genmultrand(prior.dist,prior.alpha,prior.beta,dim.draws,prior.lbnd,prior.ubnd); 

calvoPmat=zeros(dim.draws,1); 
calvoWmat=zeros(dim.draws,1); 

for jj=1:dim.draws;
    for ii=1:dim.par;
        sigmac=drawsmat(jj,pos.sigmac); 
        zeta=drawsmat(jj,pos.zeta); 
        betinDraw=drawsmat(jj,pos.betin); 
        z_ss=drawsmat(jj,pos.zss)/100; 
        bet=(betinDraw/100+1)^(-1);
        iota=drawsmat(jj,pos.iota); 
        
        calvoPmat(jj)=( 1-zeta*bet*exp( (1-sigmac)*z_ss) )*(1-zeta)/( zeta*( 1+iota*bet*exp( (1-sigmac)*z_ss ) ) ); 
        calvoWmat(jj)=( 1-zeta*bet*exp( (1-sigmac)*z_ss) )*(1-zeta)/( zeta*( 1+bet*exp( (1-sigmac)*z_ss ) ) );
                
        
    end
end

correction=(fixed.phi-1)*10 + 1; 
if flags.Kimball==1 
    calvoPmat=calvoPmat/correction; 
    calvoWmat=calvoWmat/correction;     
end 

figure; 
subplot(1,2,1); 
hist( calvoPmat,30 ); title('Slope Price PC'); 
subplot(1,2,2); 
hist( calvoWmat,30 ); title('Slope Price WC'); 


tableSlopes=tableDraws([calvoPmat calvoWmat],[],{'Slope P','Slope W'}); 
printcell(tableSlopes); 